#!/usr/bin/env python

############################################################################
#
# MODULE:	m.ecad
# AUTHOR:	Thomas Leppelt
# PURPOSE:      Import, registration,aggregation and observation of ECA&D climate data.
#
#############################################################################

#%module
#% description: Import, registration and observation of ECA&D climate data.
#% keywords: Miscellaneous
#% keywords: Climate data,ECA&D
#%end
#%option G_OPT_F_INPUT
#% description: Name of input vector file for observation
#% required : no
#% guisection: Files & format
#%end
#%option G_OPT_F_OUTPUT
#% description: Name of output vector file, merging vector map and observed ECAD climate data
#% required : no
#% guisection: Files & format
#%end
#%option
#% key: aggregation
#% type : string
#% description: Select climate data granularity for aggregation
#%options:monthly,seasonal,yearly
#% required : no
#% multiple: yes
#% answer: seasonal,yearly
#% guisection: Files & format
#%end
#%option
#% key: startcolumn
#% type : string
#% description: Name of the start date column in input file (date format: 2011-01-01)
#% required : no
#% multiple: no
#% answer: start_date
#% guisection: Files & format
#%end
#%option
#% key: endcolumn
#% type : string
#% description: Name of the end date column in input file (date format: 2011-12-31)
#% required : no
#% multiple: no
#% answer: end_date
#% guisection: Files & format
#%end
#%option
#% key: location
#% type : string
#% description: Name of location for climate analysis (compatible Lat/Lon projection required)
#% required : no
#% multiple: no
#% answer: $LOCATION_NAME
#% guisection: Files & format
#%end
#%option
#% key: observation
#% type : string
#% description: Select climate data granularity for observation
#%options:monthly,seasonal,preseasonal,yearly,preyearly,period
#% required : no
#% multiple: yes
#% answer: seasonal,yearly
#% guisection: Files & format
#%end
#%option
#% key: factor
#% type : string
#% description: Select climate factors, precipitation(rr), mean(tg)-, maximal(tx)- or minimal(tn) temperature
#%options:tg,tn,tx,rr
#% required : no
#% multiple: yes
#% answer: tg,tn,tx,rr
#% guisection: Files & format
#%end
#%flag
#% key: i
#% description: Download, import and temporal reference the ECA&D datasets into given lat/lon location
#% guisection: Projections
#%end
#%flag
#% key: o
#% description: Observe ECA&D dataset in given vector location
#% guisection: Projections
#%end

import sys
import os
import threading
import fnmatch
import grass.script as grass
from subprocess import call
import datetime
import collections

climatefactorunitconverter = {'tg' : '100.0','tn' : '100.0','tx' : '100.0','rr' : '10.0'}
calc_location = grass.gisenv()["LOCATION_NAME"]
calc_mapset = grass.gisenv()["MAPSET"]

# Download,import and temporal reference the ECA&D datasets

def ecadimport(aggregation,location,factorlist):
  
  climatefactorlist = factorlist.split(',')
  aggregationtype = {'monthly' : '1 months','seasonal' : '3 months','yearly' : '1 years'}
  aggregationmethode = {'tg' : 'average','tn' : 'average','tx' : 'average', 'rr' : 'sum'}  

  
  try:
    # Switch to given location name.
    grass.run_command('g.mapset', mapset = 'PERMANENT', location = location) 
    
  except: 
    # Create location for ECAD data processing.
    grass.run_command('g.proj', flags = 'c', epsg = 4326, location = location)
      
  for climatefactorname in climatefactorlist:
    try:
      
      # Download, unzip and import climate factors as grid ECA&D dataset in actual directory
      call(['wget', '-c', '-t', '5',str('http://eca.knmi.nl/download/ensembles/data/Grid_0.25deg_reg/'+ climatefactorname +  '_0.25deg_reg_v9.0.nc.gz')])
      call(['gunzip','-f', str(climatefactorname + '_0.25deg_reg_v9.0.nc.gz')])
      
      # Set region to 0.25 degree resolution.
      grass.run_command('g.region', nsres = 0.25, ewres = 0.25)
      
      # Import climate data
      grass.run_command('r.in.gdal', input=str(climatefactorname + '_0.25deg_reg_v9.0.nc'),
      output=str('ecad_' + climatefactorname), flags='o', overwrite = grass.overwrite)
      
      # Create txt-files of containing climate raster maps
      maplist = open(str(climatefactorname + '_map_list.txt'), "w")
      mapliststring = (grass.read_command('g.list', type='rast',flags='f').replace(' ','').split())
      del mapliststring[0:1+1]
      mapliststring = fnmatch.filter(mapliststring,str('ecad_' + climatefactorname + '*'))
      mapliststring = sorted(mapliststring,key= lambda a: float(a.split('.')[1]))
      print(mapliststring)
      for map_i in mapliststring:
          grass.run_command('r.mapcalc', expression = "%s = %s / %s" %(map_i, 
          map_i, climatefactorunitconverter[climatefactorname]), overwrite = grass.overwrite)
      mapliststring = '\n'.join(mapliststring)
      maplist.write(mapliststring)
      maplist.close()
      
      # Create and register new space time raster data sets using txt-files for daily raster sets
      grass.run_command('t.create',type='strds', output=str(climatefactorname + '_1950_2012_daily'), 
      temporal='absolute', title=str("European " + climatefactorname + " 1950-2012"),
      description=str("The european daily " + climatefactorname + " 1950 - 2012 from ECA&D"),
      overwrite = grass.overwrite)
      
      grass.run_command('t.register', flags='i', input=str(climatefactorname + '_1950_2012_daily'), 
      file=str(climatefactorname + '_map_list.txt'), start='1950-01-01', increment="1 days")
      
      #Climate data aggregation to monthly, seasonal and yearly data sets
      for aggregationlevel in aggregation.split(','):
	wherestatement = {'monthly' : str("start_time > '1950-01-01'"), 'seasonal' : str("start_time >= '1950-03-01'"), 'yearly' : str("start_time > '1950-01-01'")}
	
	grass.run_command('t.rast.aggregate', input = str(climatefactorname + '_1950_2012_daily'), 
	output = str(climatefactorname + '_1950_2012_' + aggregationlevel), 
	method = aggregationmethode[climatefactorname], 
	base = str(climatefactorname + '_' + aggregationlevel), 
	granularity = aggregationtype[aggregationlevel], overwrite = grass.overwrite, 
	where = wherestatement[aggregationlevel]) 

      #Remove downloaded file from file system 
      try:
	call(['rm', str(climatefactorname + '_0.25deg_reg_v9.0.nc')])
	
      except:
	print("# ERROR: while deleting " + str(climatefactorname) + "_0.25deg_reg_v9.0.nc #")
    except:
      print("ERROR: while processing " + str(climatefactorname) + "_0.25deg_reg_v9.0.nc.gz ")
      sys.exit(0)     

  return()
  
def ecadobserve_new(input, output, location, startcolumn, endcolumn, observation, 
                    factorlist):
    climatefactorlist = factorlist.split(',')
    aggregationmethode = {'tg' : 'average','tn' : 'average','tx' : 'average', 'rr' : 'sum'}
    
    #Switch to ECAD location
    if calc_location != location:
        grass.run_command('g.mapset', mapset = 'PERMANENT', location = location)
    oblist = observation.split(',')
    cllist = factorlist.split(',')
    
    #Split vector input by date
    inputtimelist = grass.read_command('v.db.select', map = input, 
                                        columns = str(startcolumn + ',' + endcolumn), 
                                        flags = 'c')
    inputtimelist = inputtimelist.rstrip('\n').split('\n')
    print(inputtimelist)
    for clfactor in cllist:
        print(clfactor)
        for oblevel in oblist:
            print(oblevel)
            #grass.run_command('t.vect.observe.strds', input = input, 
            #                    strds = str(clfactor + '_1950_2012_' + oblevel), 
            #                    output = str(input + '_' + clfactor + '_' + oblevel), 
            #                    vector_output = str(input + '_' + clfactor + '_' + oblevel), 
            #                    col = str(clfactor + '_' + oblevel), 
            #                    overwrite = grass.overwrite())
        count = 1
        for time in inputtimelist:
            start, end = time.split('|')
            start_str = str(start + ' 00:00:00')
            end_str = str(end + ' 00:00:00')
            time_str = str("start_time = '" + start_str + "'")
            print(time_str)
            grass.run_command('t.vect.db.select', col = str(clfactor + '_' + oblevel),
                                input = str(input + '_' + clfactor + '_' + oblevel), 
                                t_where = time_str, where = str("cat = " + str(count)))
            count = count + 1
                                
    return()
    
def ecadobserve(input, output, location, startcolumn = "start_date", 
                endcolumn = "end_date", observation = "seasonal,yearly",
                factorlist = "tg,tn,tx,rr"):
  
  climatefactorlist = factorlist.split(',')
  aggregationmethode = {'tg' : 'average','tn' : 'average','tx' : 'average', 'rr' : 'sum'}
  
  #Switch to ECAD location
  grass.run_command('g.mapset', mapset = 'PERMANENT', location = location)
  
  #Create space time vector dataset for vector input
  grass.run_command('t.create', type = 'stvds', output = input, title = 'Vector input for climate analysis', 
  description = 'Vector input for climate analysis',overwrite = grass.overwrite)
  
  #Split vector input by date
  inputstartlist = grass.read_command('v.db.select', map = input, columns = startcolumn, flags = 'c')
  inputstartlist = inputstartlist.rstrip('\n').split('\n')
  
  inputendlist = grass.read_command('v.db.select', map = input, columns = endcolumn, flags = 'c')
  inputendlist = inputendlist.rstrip('\n').split('\n')
  
  inputmergelist = zip(inputstartlist,inputendlist)
  inputmergelist = collections.OrderedDict.fromkeys(inputmergelist).keys()
  inputfilelist =[]
  nextract=0
  for starti,endi in inputmergelist:
    grass.run_command('v.extract', input = input, 
    where = str(startcolumn + " = '" + starti + "' AND " + endcolumn + " = '" + endi + "'"),
    output = str(input + "_" + str(nextract)), overwrite = grass.overwrite)
    inputfilelist.append(str(input + "_" + str(nextract)))
    nextract+=1
  
  inputmaplist = open(str(input + '_map_list.txt'), "w")

  for iextract in range(0,nextract):
    registerinputline = str(inputfilelist[iextract] +'|' + '|'.join(inputmergelist[iextract]) + '\n')
    inputmaplist.writelines(registerinputline)
  inputmaplist.close()
  
  #Register vector data in space time vector dataset
  grass.run_command('t.register', input = input, type = 'vect', file = str(input + '_map_list.txt'), 
  sep = '|', flags='i')

  #create maplist for previous year
  formatlist = []  
  formatlist[1:len(inputstartlist)] = ["%Y-%m-%d"] * len(inputstartlist)
  inputstartdatelist = map(datetime.datetime.strptime, inputstartlist, formatlist)
  for dateelement in range(0,len(inputstartdatelist)):
    inputstartdatelist[dateelement] = inputstartdatelist[dateelement].replace(
			      year = inputstartdatelist[dateelement].year - 1)
    inputstartdatelist[dateelement] = datetime.datetime.strftime(inputstartdatelist[dateelement],"%Y-%m-%d")
  
  formatlist = []  
  formatlist[1:len(inputendlist)] = ["%Y-%m-%d"] * len(inputendlist)
  inputenddatelist = map(datetime.datetime.strptime, inputendlist, formatlist)
  for dateelement in range(0,len(inputenddatelist)):
    inputenddatelist[dateelement] = inputenddatelist[dateelement].replace(
			      year = inputenddatelist[dateelement].year - 1)
    inputenddatelist[dateelement] = datetime.datetime.strftime(inputenddatelist[dateelement],"%Y-%m-%d")

  grass.run_command('t.create', type = 'stvds', output = str(input + '_previous'), title = 'Previous year vector input for climate analysis', 
  description = 'Previous year vector input for climate analysis',overwrite = grass.overwrite)
  
  inputmergedatelist = zip(inputstartdatelist,inputenddatelist)
  inputmergedatelist = collections.OrderedDict.fromkeys(inputmergedatelist).keys()
  inputfiledatelist =[]
  nextract=0
  for starti,endi in inputmergelist:
    grass.run_command('v.extract', input = input, 
    where = str(startcolumn + " = '" + starti + "' AND " + endcolumn + " = '" + endi + "'"),
    output = str(input + "_previous_" + str(nextract)), overwrite = grass.overwrite)
    inputfiledatelist.append(str(input + "_previous_" + str(nextract)))
    nextract+=1
    
  inputmapdatelist = open(str(input + '_map_list_previous.txt'), "w")

  for iextract in range(0,nextract):
    registerinputdateline = str(inputfiledatelist[iextract] +'|' + '|'.join(inputmergedatelist[iextract]) + '\n')
    inputmapdatelist.writelines(registerinputdateline)
  inputmapdatelist.close()
  
  #Register vector data in space time vector dataset for previous year
  grass.run_command('t.register', input = str(input + '_previous'), type = 'vect', file = str(input + '_map_list_previous.txt'), 
  sep = '|', flags='i')
    
  for observationlevel in observation.split(','):
    
    for climatefactorname in climatefactorlist:   
      #Observation of climatedata
      if observationlevel == 'seasonal':
	seasonstartmonth = {'spring' : "'2 months'",'summer' : "'5 months'",'autumn' : "'8 months'",'winter' : "'11 months'"}
	
	for seasonname in ['spring','summer','autumn','winter']:	  
	  grass.run_command('t.rast.extract', overwrite = grass.overwrite, 
	  input = str(climatefactorname + '_1950_2012_' + observationlevel), 
	  output = str(climatefactorname + '_1950_2012_' + seasonname), 
	  where = str("start_time = datetime(start_time,'start of year', " + seasonstartmonth[seasonname] + ")"),
	  #expression = str(climatefactorname + '_1950_2012_' + observationlevel + '/' + climatefactorunitconverter[climatefactorname]), 
	  base = str(climatefactorname + '_1950_2012_' + seasonname))
	  
	  grass.run_command('t.vect.what.strds', input = input, 
	  strds = str(climatefactorname + '_1950_2012_' + seasonname),
	  overwrite = grass.overwrite, column = str(climatefactorname + '_' + seasonname),
	  sampling='during,overlap,contain,equal')
	  
	  for iextract in range(0,nextract):
	    columnlist = grass.vector_columns(str(input + "_" + str(iextract))).keys()
	    climatecolumnname = str(climatefactorname + '_' + seasonname)
	    if climatecolumnname not in columnlist:
	      grass.run_command('v.db.addcolumn', map = str(input + "_" + str(iextract)), 
	      columns = str(climatefactorname + '_' + seasonname + ' double precision'))
	      grass.core.warning(str(input + "_" + str(iextract) + ' missing value for ' + climatefactorname + '_' + seasonname))
      
      if observationlevel == 'preseasonal':	
	seasonstartmonth = {'spring' : "'2 months'",'summer' : "'5 months'",'autumn' : "'8 months'",'winter' : "'11 months'"}
	
	for seasonname in ['spring','summer','autumn','winter']:	  
	  grass.run_command('t.rast.extract', overwrite = grass.overwrite, 
	  input = str(climatefactorname + '_1950_2012_' + 'seasonal'), 
	  output = str(climatefactorname + '_1950_2012_' + seasonname), 
	  where = str("start_time = datetime(start_time,'start of year', " + seasonstartmonth[seasonname] + ")"),
	  #expression = str(climatefactorname + '_1950_2012_' + 'seasonal' + '/' + climatefactorunitconverter[climatefactorname]), 
	  base = str(climatefactorname + '_1950_2012_' + seasonname))
	  
	  grass.run_command('t.vect.what.strds', input = str(input + '_previous'), 
	  strds = str(climatefactorname + '_1950_2012_' + seasonname),
	  overwrite = grass.overwrite, column = str(climatefactorname + '_pre' + seasonname),
	  sampling='during,overlap,contain,equal')
	  
	  for iextract in range(0,nextract):
	    columnlist = grass.vector_columns(str(input + "_previous_" + str(iextract))).keys()
	    climatecolumnname = str(climatefactorname + '_pre' + seasonname)
	    if climatecolumnname not in columnlist:
	      grass.run_command('v.db.addcolumn', map = str(input + "_previous_" + str(iextract)), 
	      columns = str(climatefactorname + '_pre' + seasonname + ' double precision'))
	      grass.core.warning(str(input + "_" + str(iextract) + ' missing value for ' + climatefactorname + '_pre' + seasonname))

	    grass.run_command('v.db.addcolumn', map = str(input + "_" + str(iextract)), 
	    columns = str(climatecolumnname + ' double precision'))
	  
	    grass.run_command('db.execute', sql = str('UPDATE ' + input + "_" + str(iextract) + 
	    ' SET ' + climatecolumnname + '=(select ' + climatecolumnname + ' from ' + input + 
	    "_previous_" + str(iextract) + ' where ' + input + "_" + str(iextract) + '.cat=' + 
	    input + "_previous_" + str(iextract) + '.cat)')) 
	  
      if observationlevel == 'yearly':	
	grass.run_command('t.rast.extract', overwrite = grass.overwrite, 
	input = str(climatefactorname + '_1950_2012_' + observationlevel), 
	output = str(climatefactorname + '_1950_2012_' + 'annual'), 
	#expression = str(climatefactorname + '_1950_2012_' + observationlevel + '/' + climatefactorunitconverter[climatefactorname]), 
	base = str(climatefactorname + '_1950_2012_' + 'annual'))
	
	grass.run_command('t.vect.what.strds', input = input, 
	strds = str(climatefactorname + '_1950_2012_' + 'annual'),
	overwrite = grass.overwrite, column = str(climatefactorname + '_annual'),
	sampling='during,overlap,contain,equal')
	
	for iextract in range(0,nextract):
	  columnlist = grass.vector_columns(str(input + "_" + str(iextract))).keys()
	  climatecolumnname = str(climatefactorname + '_annual')
	  if climatecolumnname not in columnlist:
	    grass.run_command('v.db.addcolumn', map = str(input + "_" + str(iextract)), 
	    columns = str(climatefactorname + '_annual' + ' double precision'))
	    grass.core.warning(str(input + "_" + str(iextract) + ' missing value for ' + climatefactorname + '_annual'))
	    
      if observationlevel == 'preyearly':	
	grass.run_command('t.rast.extract', overwrite = grass.overwrite, 
	input = str(climatefactorname + '_1950_2012_' + 'yearly'), 
	output = str(climatefactorname + '_1950_2012_' + 'annual'), 
	#expression = str(climatefactorname + '_1950_2012_' + 'yearly' + '/' + climatefactorunitconverter[climatefactorname]), 
	base = str(climatefactorname + '_1950_2012_' + 'annual'))
	
	grass.run_command('t.vect.what.strds', input = str(input + '_previous'), 
	strds = str(climatefactorname + '_1950_2012_' + 'annual'),
	overwrite = grass.overwrite, column = str(climatefactorname + '_preannual'),
	sampling='during,overlap,contain,equal')
	
	for iextract in range(0,nextract):
	  columnlist = grass.vector_columns(str(input + "_previous_" + str(iextract))).keys()
	  climatecolumnname = str(climatefactorname + '_preannual')
	  if climatecolumnname not in columnlist:
	    grass.run_command('v.db.addcolumn', map = str(input + "_previous_" + str(iextract)), 
	    columns = str(climatefactorname + '_preannual' + ' double precision'))
	    
	    grass.core.warning(str(input + "_previous_" + str(iextract) + ' missing value for ' + climatefactorname + '_preannual'))
	  
	  grass.run_command('v.db.addcolumn', map = str(input + "_" + str(iextract)), 
	  columns = str(climatecolumnname + ' double precision'))
	  
	  grass.run_command('db.execute', sql = str('UPDATE ' + input + "_" + str(iextract) + 
	  ' SET ' + climatecolumnname + '=(select ' + climatecolumnname + ' from ' + input + 
	  "_previous_" + str(iextract) + ' where ' + input + "_" + str(iextract) + '.cat=' + 
	  input + "_previous_" + str(iextract) + '.cat)')) 
	  
      if observationlevel == 'monthly':	
	monthnamelist = {'january' : "'0 months'",'february' : "'1 months'",'march' : "'2 months'",'april' : "'3 months'",'may' : "'4 months'",'june' : "'5 months'",'july' : "'6 months'",'august' : "'7 months'",'september' : "'8 months'",'october' : "'9 months'",'november' : "'10 months'",'december' : "'11 months'"}
	
	for monthname in ['january','february','march','april','may','june','july','august','september','october','november','december']:	  
	  grass.run_command('t.rast.extract', overwrite = grass.overwrite, 
	  input = str(climatefactorname + '_1950_2012_' + observationlevel), 
	  output = str(climatefactorname + '_1950_2012_' + monthname), 
	  where = str("start_time = datetime(start_time,'start of year', " + monthnamelist[monthname] + ")"),
	 # expression = str(climatefactorname + '_1950_2012_' + observationlevel + '/' + climatefactorunitconverter[climatefactorname]), 
	  base = str(climatefactorname + '_1950_2012_' + monthname))
	  
	  grass.run_command('t.vect.what.strds', input = input, 
	  strds = str(climatefactorname + '_1950_2012_' + monthname),
	  overwrite = grass.overwrite, column = str(climatefactorname + '_' + monthname),
	  sampling='during,overlap,contain,equal')
	  
	  for iextract in range(0,nextract):
	    columnlist = grass.vector_columns(str(input + "_" + str(iextract))).keys()
	    climatecolumnname = str(climatefactorname + '_' + monthname)
	    if climatecolumnname not in columnlist:
	      grass.run_command('v.db.addcolumn', map = str(input + "_" + str(iextract)), 
	      columns = str(climatefactorname + '_' + monthname + ' double precision'))
	      grass.core.warning(str(input + "_" + str(iextract) + ' missing value for ' + climatefactorname + '_' + monthname))

      if observationlevel == 'period':	
	grass.run_command('t.rast.aggregate.ds', input = str(climatefactorname + '_1950_2012_daily'), 
	output = str(climatefactorname + '_1950_2012_' + observationlevel),
	sample = input, type = 'stvds',
	method = aggregationmethode[climatefactorname], 
	base = str(climatefactorname + '_' + observationlevel), overwrite = grass.overwrite) 

	grass.run_command('t.rast.extract', overwrite = grass.overwrite, 
	input = str(climatefactorname + '_1950_2012_' + observationlevel), 
	output = str(climatefactorname + '_1950_2012_' + 'periodical'), 
	#expression = str(climatefactorname + '_1950_2012_' + observationlevel + '/' + climatefactorunitconverter[climatefactorname]), 
	base = str(climatefactorname + '_1950_2012_' + 'periodical'))
	
	grass.run_command('t.vect.what.strds', input = input, 
	strds = str(climatefactorname + '_1950_2012_' + 'periodical'),
	overwrite = grass.overwrite, column = str(climatefactorname + '_period'))
  
  if len(inputfilelist) != 1:
    grass.run_command('v.patch', input = ','.join(inputfilelist), output = output, 
    flags = 'e', overwrite = grass.overwrite)
  else:
    grass.run_command('v.extract', input = inputfilelist, output = output, overwrite = grass.overwrite)
    
  return()

def main():
    input = options['input']
    output = options['output']
    aggregation = options['aggregation']
    startcolumn = options['startcolumn']
    endcolumn = options['endcolumn']
    location = options['location']
    observation = options['observation']
    factorlist = options['factor']
    ecadimpflag = flags['i']
    ecadobservflag = flags['o']
    
    if location == "$LOCATION_NAME":
        envdict = grass.gisenv()
        location = envdict["LOCATION_NAME"]
        
    if ecadimpflag:
      ecadimport(aggregation,location,factorlist)
      
    if ecadobservflag:
      ecadobserve(input,output,location,startcolumn,endcolumn,observation,factorlist)

if __name__ == "__main__":
    options, flags = grass.parser()
    main()
